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The quantum transport through nanoscale junctions is governed by the charging energy U of 
the device. We employ the recently developed scattering-states numerical renormalization group 
approach to open quantum systems to study nonequilibrium Green's functions and current- voltage 
characteristics of such junctions for small and intermediate values of U. We establish the accuracy 
of the approach by a comparison with diagrammatic Kadanoff-Baym-Keldysh results which become 
exact in the weak coupling limit U — > 0. We demonstrate the limits of the diagrammatic expansions 
at intermediate values of the charging energy. While the numerical renormalization group approach 
correctly predicts only one single, universal low-energy scale at zero bias voltage, some diagrammatic 
expansions yield two different low-energy scales for the magnetic and the charge fluctuations. At 
large voltages, however, the self-consistent second Born as well as the GW approximation reproduce 
the scattering-states renormalization group spectral functions for symmetric junctions, while for 
asymmetric junctions the voltage-dependent redistribution of spectral weight differs significantly 
in the different approaches. The second-order perturbation theory does not capture the correct 
single-particle dynamics at large bias and violates current conservation for asymmetric junctions. 

PACS numbers: 73.21.La, 73.63. Rt, 72.15.Qm 



I. INTRODUCTION 

Quantum dots and single-molecule junctions have 
been considered as possible building blocks for nano- 
electronics and quantum-information processingP^ R e _ 
cent technological progress has made it possible to man- 
ufacture and study electron transport trough ultra-small 
quantum-dot devices, nanotubes or single moleculesP^J] 
These devices are designed as a small central region com- 
prising of a quantum dot or a single molecule which is 
coupled to at least two leads where the finite bias volt- 
age is applied to. The investigation of such devices is of 
fundamental importance for our understanding of open 
quantum systems out of equilibrium. 

Due to the quantization of the charge, the physical 
properties of such junctions are dominated by many- 
body effects at temperatures below the charging energy 
U = e 2 /(2C), where C is the capacitance of the device. 
The experimental devices are often fully controllable by 
external gate electrodes or elongation of the scanning 
tunneling microscope tipP This gives the opportunity to 
directly study true many-body correlation effects, such as 
the Kondo effect (see for example Ref . and [T3|) . under 
the influence of an external bias voltage. However, the 
theoretical understanding of the interplay between co- 
herent transport favored by many-body correlations and 
current-driven dephasing at finite bias is still at its in- 
fancy and further investigations are needed. 

The present work has two main objectives. On the one 
hand, we establish the reliability of the recently intro- 
duced scattering-states numerical renormalization group 
(SNRG) approactP^ to quantum transport by compar- 
ing results for small values of U to the diagrammatic 
Kadanoff-Baym-Keldysh approach, which becomes exact 



in the limit U — > 0. On the other hand, we will discuss 
discrepancies and reveal shortcomings of those diagram- 
matic approaches at intermediate values of the charging 
energy. 

We investigate quantum transport through a quantum- 
dot device using a minimal modeP^ where the complex 
interacting region is replaced by a single spinful orbital 
which is coupled to two noninteracting leads. A single 
Coulomb matrix element U accounts for the charging en- 
ergy of the device. We calculate nonequilibrium spectral 
functional and current-voltage (IV) characteristics using 
the SNRG as well as different approximations^"^ within 
the diagrammatic Kadanoff-Baym-Keldysh expansion in 
the local Coulomb interaction U. 

Over the past 40 years, the Keldysh technique^ has 
proven to be the most successful approach to nonequi- 
librium dynamics. In the context of quantum trans- 
port through nano-junctions direct expansions in the 
interactiorP2HHl as we ll as self-consistent re-summation 
schemes have been employed Q3ES1 However, such dia- 
grammatic expansions rely on a small expansion param- 
eter, and are, therefore, confined to weak coupling. But 
quantum-impurity models^ commonly used in the the- 
ory of quantum transport on the molecular level often 
exhibit infra-red divergences in perturbation theory^ 
which also restrict the diagrammatic Keldysh approaches 
to certain parameter regimes usually to high temperature 
or to large bias. 

In contrast to equilibrium conditions, where complete 
and accurate solutions can be obtained using a vari ety of 
nonperturbative techn iques such as the Bethe ansatzf^M 
conformal field theorypSM or Wilso n's n umerical renor- 
malization group (NRG) approach] 26 * 31 ! techniques for 
calculating quantum-transport out of equilibrium re- 
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main largely at the development stage. Recent advance- 
ments on the analytical sid e 32 " 37 ! incl ude s uitable adap- 
tations of the Wegner's flow-equation 38 39 and the real- 
time renormalization-group methodJ^EH These methods 
can successfully access large voltages, but are generally 
confined to the weak-coupling regime. Base d on the 
scattering-states approach to quantum transport! 36 ! 43 ! the 
Bethe ansatz was extended to quantum-impurity models 
out of equilibrium,^ but remains limited to a certain 
class of models. 

On the numerical side, progress has been made in sev- 
eral directions. Currents have been extracted from time- 
dependent density matrix renormalization group^^ cal- 
culations using finite ID wires, and the results agree well 
with Bethe ansatz results for certain modelsP^ Quantum 
Monte Carlo approaches based on scat tering states^ can 
access the intermediate coupling regime^ESl at finite bias. 
Recent real-time formulations of continuous-time quan- 
tum Monte CarlcES^l and an iterative real-time path in- 
tegral approach^ to quantum transport offer the appeal- 
ing advantage of working directly in the continuum limit, 
but are confined to relatively short time scales. Access 
to low temperatures and long times is hampered in the 
former case by a severe sign problem, and by the extrapo- 
lation to long memory times in the latter case. Hence nei- 
ther approach can presently be applied to nonequilibrium 
dynamics of correlated systems with a small underlying 
energy scale, as is the case with ultra-small quantum dots 
when tuned to the Kondo regime. 

The usage of Lippmann-Schwinger scattering states 
has been well established in quantum-field theory^ for 
over 50 years and also successfully adapted to the descrip- 
tion of quantum transport through stro ngly inter acting 
nano-devices coupled to ballistic leads |M3 7 | 43 | 56| These 
states fulfill the correct boundary condition of the open 
quantum system: (i) they break time-reversal symmetry 
and, therefore, are (ii) complex and current-carrying and 
(iii) describe ballistic transport in the leads combined 
with scattering events in the small interacting quantum- 
dot region. This time-reversal symmetry breaking is re- 
quired for current carrying systems and reveals itself nat- 
urally in all diagrammatic approaches by the occurrence 
of retarded and advanced Green's functions. It is a con- 
sequence of any regularization when performing the limit 
to an infinitely large system. 

In particular, the work of HershfielcP^ and Doyon 
and AndreP^ has rigorously shown that these boundary 
conditions remain unaltered when a local interaction is 
switched on. The noninteracting current-carrying sys- 
tem evolves into the new steady-state of the interacting 
system, and the steady-state density operator retains a 
Boltzmannian formP^l The explicit construction of those 
scattering states allows to exactly solve the DC and AC 
Kondo model at the Toulouse pohrP ^ED as we ll as the 
interacting resonant level model.^EU 

Recently, an extension^ to Wilson's numerical renor- 
malization group has been developed for steady-state 
quantum transport through nano-devices which is able 



to deal with the crossover from weak to strong coupling 
for arbitrary bias voltages. It is based on Oguri's idea 60 
of discretizing the single-particle scattering states which 
are the solutions of the Lippman-Schwinger equation^ 
for the noninteracting problem and, therefore, fulfill the 
correct boundary condition of an open quantum system. 
This scattering-states numerical renormalization group 
approactP31 (SNRG) evolves the analytically known den- 
sity matrix of a noninteracting system to the density ma- 
trix of the fully interacting prob lem b y employing the 
time-dependent NRG (TD-NRG).EMI] The NRG is ide- 
ally suited to the problem, being known to provide accu- 
rate solutions of quantum- impurity models on all relevant 
interaction strengths at zero biasP^ Sin ce th e TD-NRG 
can access exponentially long time scales EMU dwell times 
on the order of the inverse Kondo-temperature are easily 
accessible. 

This paper is organized as follows. After the model 
used is defined, we provide the details of the different 
theoretical approaches in Sec. [IT] We summarize the ba- 
sic ideas of the SNRG method introduced in Ref. [14] in 
Sec. |II B| and state all necessary equations of the dia- 
grammatic nonequilibrium techniques in Sec. |II C| The 
main body of the paper is in Sec. Ill where we present 



and discuss the results obtained for the various meth- 
ods. In order to set the stage for a detailed comparison 
between the SNRG and diagrammatic approaches at fi- 
nite bias, we begin with a discussion of the magnetic and 
charge fluctuation scales at zero bias in Sec. |III A| Since 
the NRG provides an accurate solution in this regime for 
arbitrary coupling strengths and temperatures, this re- 
veals the validity range of the diagrammatic approaches. 
We show that — in contrast to the NRG — some of the 
diagrammatic expansions fail to produce a single low- 
energy scale for intermediate and large values of U . How- 
ever, in the weak correlation regime all these approaches 
agree excellently for arbitrary voltages at small U and 
yield the same nonequilibrium Green functions as well 
as IV characteristics which are presented in Sec. IIIB 



Discrepancies between the different approaches at inter- 
mediate values of the Coulomb interaction are discussed 
in Sec. |III C[ where the spectral functions and IV curves 
of a symmetric and an asymmetric junction are consid- 
ered. We conclude with summary and a short outlook in 
Sec. HVl 



II. THEORY 
A. Model 

Quantum impurity models are used to describe quan- 
tum transport on the molecular level. Their Hamiltonian 
U 



rl — Wimp + Wbath + Hi 



(1) 



consists of three parts: an impurity part T-Li mv modeling 
the interacting device with a finite number of degrees of 
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freedom, one or several bosonic or fermionic baths rep- 
resented by Hbath, and the coupling of these subsystems 

Throughout this paper, we restrict ourselves to junc- 
tions modeled by the single impurity Anderson model 
with one spinful orbital coupled to a left (L) and a right 
(R) lead and an on-site Coulomb repulsion U 
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Here, Ed is the single-particle energy of the quantum 
dot, h d a — S a d a measures its orbital occupancy and t acr 
represent the elementary hybridization-matrix elements 
coupling the dot to the two leads. The different chemical 
potentials p a in both leads appear as a shift of the band 
centers and are functions of the external voltage V — 
PR - PL- 

For simplicity, we assume that both leads have the 
same density of states, p R (e) = Pl( 6 ) = p( e ); charac- 
terized by the same band width D but different band 
centers. This Hamiltonian is commonly used to model a 
single Co ulom b-blockade resonance in ultra-small quan- 
tum dots.EESl 



B. Scattering-states numerical renormalization 
group approach 

1. Definition of the scattering states 

In the absence of the local Coulomb repulsion Hu — 
Un^n^, the single-particle problem is diagonalized ex- 
actly in the continuum Urmll 14 ' 43 ' 4iJ ' 56 ' tiu ' B3 ' 64 ' by the fol- 
lowing scattering-states creation operators 



7La = c L a + t a \/pJejG r 0lT {e) 



\p f , J~t<x'\/P<x'{z') t 

/ e + iS- e ' t'™' 



(3) 



a = L(R) labels left (right) moving scattering states 
The local retarded resonant level 



created by ll aL[R) - 
Green's function 
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enters as an expansion coefficient. Defining t 
\Jt\ + t 2 R , we will use rR(£) = t R ^/t and 



AH = ^rljde 



Pa(e) 
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in the following. 

In the limit of infinitely large leads — volume Vol. — > 
oo — the single-particle spectrum remains unaltered, and 
these scattering states diagonalize the Hamiltonian ^ 
for U = 0: 



Ht = U{U = 0) = 53 f dee 7 l a7e<7Ct 

T TD. _ ** 
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The scattering states are solutions of the Lippmann- 
Schwinger equation^ and therefore break time-reversal 
symmetry, which constitutes a necessary boundary con- 
dition to describe a current carrying open quantum sys- 
tem. This is encoded in the small imaginary part +id~ 
entering Eq. ([3]) - ^ required for convergence when per- 
forming the continuum limit Vol. — > oo in the leads. 

The complex expansion coefficients in ^ are given by 
retarded functions, e.g. Gg CT (e), which causes the scatter- 
ing states to be complex and current carrying. For zero 
bias voltage, time-reversal symmetry manifests itself in 
the identical spectrum for left and right movers which 
are time-reversal pairs in that limit. 

To avoid any contribution from bound states, we 
will implicitly assume a wide band limit: D 3> 
max{\E d \,r, \V\}, where Y a = nt 2 a p{0) and T = T L +T R . 

Hershfield has shown that the density operator for such 
a noninteracting current carrying quantum system re- 
tains its Boltzmannian forrrP^ 

PO = f „,„„. ,-. ,1 > ^0 = V Pa / de 7Lc7e<raC0 
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even for finite bias. The Yq operator accounts for the dif- 
ferent occupation of the left- and right-moving scattering 
states, and p a for the different chemical potentials of the 
leads. 

Therefore, all steady-state expectation values of oper- 
ators can be calculated using po which includes the finite 
bias. In the absence of a Coulomb repulsion U , this is a 
trivial and well-understood problem. It was shown 60 that 
the current expectation value usi ng this d ensity-operator 
po reproduced the standard resuh r 5 ! 65 ! 66 ! for noninteract- 
ing devices. The knowledge of the analytical form of po, 
however, make s this steady-state model accessible to a 
NRG approachPHU 

The expansion coefficients of ^\ aa in Eq. ^ contain 
the complex single-particle Green function G , Q CT (e) which 
we separate in modulus and phase 



G r At) = |GUe)|e-^) 



(8) 



This phase is absorbed into the new scattering states 
itaa -> itaa = iLa^"^ by a local gauge transforma- 
tion. The impurity operator d\ is expanded into left- and 
right-mover contributions 



4 = r R dl R + r L dl 



(9) 
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using the inversion of Eq. ([3|. These two new operators 
dl„ are then defined as 



4 a = i / deVpM|G5 CT (e)|7L 



(10) 



and obey the anti-commutator relation {d aa >d a , _,} 



cretization of the scattering-states continuum 7 £(JQ in 
intervals J£ = [A-^+^D, A^^+^^Z?] and 7™ = 
[_A-(n+ z -i)) j D ! _ A-(" +^)D1 ( n = 1,2, •••), controlled 
by the parameter d 26 ' 311 A > 1 and z £ (0,1]. The in- 
tervals for n — are defined as 19 = [A~ 2 D,.D] and 
7° = [— D, — A~ Z D]. An average over various z-values^ 
is used to mimic the conduction band continuum. 



2. Discretization of the scattering states 

The scattering-states numerical renormalization group 
approach^ (SNRG) starts from a logarithmic dis- 



Then, the discretized version of the noninteracting 
Hamiltonian ([6| is mapped onto a semi-infinite Wilson 
chain 



oo oo 
Ho(A) = y ]y ] Wncrafnaafna-a + ^ ] (tncrafnaafn+lcra + ^ncrafn+laafnua) (H) 
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whose tight-binding matrix elements t ntTa decay expo- 
nentially t naa ex K~ n / 2 for large n. In contrast to the 
standard NRGf^ES] the impurity degree of freedom has 
been included into Hq(A) since not the leads but the 
full scattering states have been discretized. Any com- 
plex phase in the tight-binding parameters t naa can be 
absorbed into the creation (anihilation) operators f\ aa 
ifnaa) of an electron on the chain link n with spin a and 
mover a by a local gauge transformation. 

We use d aa defined in Eq. (10) as starting vector 
foaa — dcra for the Householder transformation^ and 
obtain the tight-binding coe fficien ts of the Wilson chain 
(111 by the usual procedure! 26 * 31 ! It is straight forward 



to shown that the energy of the first chain link corre- 
sponds to the energy of the original quantum-dot orbital: 



WQa 



E ri . 



3. Local Coulomb interaction 

In order to include the local Coulomb interaction, the 
density operator n d = d\d a must be expanded in the new 
orbitals d aa . It consist of two contributions: A density 



term and a backscattering term h d a 



O 



back 
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and the backscattering 0^ ack term is defined as 

6 b a ack = r L r R (dl R d aL + di L d aR ) . (13) 
The local Coulomb interaction Hfj 

H v = U ^nJ + ^O^ acfc ^ CT + 0| acfc df c ^(14) 



leads to a mixing of left and right movers since 0^ ack 
does not commute with Yq. However, the term Hy, 



U 



- 1 



(15) 



commutes with Y Q and can be absorbed into the steady- 
state density operator p n — > p = exp[— f3{W — Y )]/Z 
with W — Hq + Hjj using the arguments given in Ref.l36l 



4- Review of the time- dependent numerical renormalization 
group approach 

Starting from an equilibrated system for times t < 0, 
the initial Hamiltonian % is changed to 77/ by a sudden 
quench at t — 0. Then, the density operator p(t) evolves 
from its initial value po at t — as 



e -i-H f t/n^ e m f t/h 



(16) 



If tti(f) describes a quantum impurity problem and O is 
an impurity operator, it was recently shown that the real- 
time dynamics o f the expectation value of 0(t) = (0(t)) 
can be calculate d 61 1 62 1 by evaluating 



(lis 



0(t) 



VV 



red 

/ , / , Pr,s 



(m)o; 
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where 0™ r = (s,e;m\0\r,e;m) denotes the matrix ele- 
ments of the operator O and E™ the NRG eigenenergy of 
the eigenstate |r; m) to 77/ at NRG iteration m. The sum 



restriction Ylrs indicates that at least one of the states 
r, s must be a discard state at iteration m. Excitations 
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between two retained states will be refined in the follow- 
ing iterations m! > m and, therefore, will contribute at a 
later iteration, e labels the environment degrees of free- 
dom of the Wilson chain links to be incorporated in sub- 
sequent iterations m! > m and the \r, e; m) = |r; m) ® \e) 
are just tensor-product states of the eigenstates of the 
mth iteration and the yet uncoupled rest chain. The re- 
duced density matrix 



red 
rr.s 



(rn) 



E< 



r, e; m\po\s, e; m) 



(18) 



traces out all environment degrees of freedom e. The ini- 
tial conditions are encoded into the density operator po 
calculated with the initial Hamiltonian Hi. The calcula- 
tion of the overlap matrix between the NRG eigenstates 
of Hi and Hf allows for the basis set transformation of 
p r r e f{m) into the basis of the final Hamiltonian provided 
that po remains restricted to the last Wilson shell jSUSS 
This transformed p r r e f(m) enters Eq. (17). 

The discarded states form a complete basis selP^ for 
the Fock-space of the entire Wilson chain of length N, i.e. 
Fm = span{|Z, e; m}} where I labels all discarded states at 
iteration m. The iterative diagonalization thus procures 
the set of (approximate) eigenstates for the whole energy 
range from high energies on the order of the bandwidth 
down to very low energies such as the Kondo scale. This 
is indispensable because nonequilibrium processes usu- 
ally involve all energy scales and cannot be confined to a 
finite low energy window set by the last Wilson shell as 
in the usual equilibrium NRG. 



The scattering-states NRG approach and steady-state 
Green's function 



In Sec. II B 1 we have argued that the analytic form the 
steady-state nonequilibrium density operator is known 
for the noninteracting case. This allows for applying the 
NRG approach to construct a faithful representation of 
po (V, U = 0) . We assume that when switching on the 
Coulomb interaction Hjj for infinitely large leads (i) a 
steady state is reached after some characteristic but finite 
time and (ii) it is unique and independent of the initial 
condition. As described earlier, the boundary condition 
of time-reversal symmetry breaking is imposed on the 
scattering states and the nonequilibrium density operator 
at t = for U — 0. The interaction quench at t = 0, i.e. 
switching on a local scattering potential, and the subse- 
quent unitary time evolution do not affect this boundary 
condition, and the time-evolved operators characterize 
the interacting current carrying open quantum system. 

The time average of the density operator 



Po 



= lim 



1 



dtp(t) 



(19) 



projects out the steady-state contributions 
to the time-evolved density operator p(t) — 



cxp(—iHft/h)pQexp(iHft/h) even in a finite-size 

system: only the energy diagonal terms contribute in 

accordance with the steady-state condition [Hf, Poo] = 0. 

Even though p^ remains unknown analytically, we can 

cons truct it systematically using the time-dependent 
NRG f6il£2] described above 

The steady-state retarded Green's function is defined 

as 



G\ B {t) = -iTr 



p 00 [A{t),B] s e(t) 



(20) 



where A(t) = e m f t f n Ae- m t t l h , [A(t),B] s denotes the 
commutator (s = —1) for bosonic, and the anti- 
commutator (s = 1) for fermionic correlation functions. 
This Green's fu nction can be calculated using the time- 
dependent NR d 61 | 62 J and extending ideas developed for 
equilibrium Green's functions.^ The completeness rela- 
tion for the basis of discarded states as introduced above 
is given by 



N 



1 = E EEi^ m M 



e; m\ 



(21) 
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where m m i n denotes the first iteration at which the NRG 
truncation is employed, N is the total number of itera- 
tions (i.e. the length of the Wilson chain), and I only 
runs over the states which are discarded at iteration m. 
For each ite ration m, we can partition the completeness 
relation ( 21 ) into two parts, 1 = l m + l+, where the first 
part incorporates the iterations m m i n to m and the sec- 
ond the iterations m + 1 to N. Since 1+ spans the part 
of the Fock-space which contains all kept states \k,e;m) 
after iteration m, the identity 

N 

l m = E E E \l,e;m'){l,e;m'\ 

m'— m+1 l£dis e 
k£kept e 

must hold. The different contributions to the Green's 
function are calculated for each energy scale D m oc 
A~ m / 2 by expanding the (anti-)commutator in Eq. (20) 



and inserting the completeness relations Eqs. (21) and 
( 22 ) repeatedly. By making use of the fact that local 



operators A and B are diagonal in the environment de- 
gree of freedom e, reduced density matrices p^ e f(m) oc- 
cur naturally when tracing out the environment e here as 
well. Although the excitation energies remain confined to 
the same energy scale, terms connecting different energy 
scales D m and D m i are implicitly included through the 
reduced density matrices such as defined in (18). Simi- 



lar to the real-time dynamics, the summation over all m 
then ensures that all en ergy scales D m contribute to the 
Green's functions! 16 * 68 ! A detailed derivation is given in 
Ref. [TSJ It was shown that the algorithm is identical to 
the equilibrium algorithrrPO if Hi = Hf. Laplace trans- 
forming G r A B (t) yields the steady-state spectral function 
for the retarded Green's function which is used to calcu- 



late the current (see Eq. ( 43 ) below) 



G 



C. The Kadanoff-Baym-Keldysh approach 

We employ the nonequilibrium perturbation theory 
as formulated by Kadanoff and BayrrP*] and KeldystPH 
on the usual Keldysh time contour, for example, see 
Refs. I70HT21 Since we are only interested in the steady- 
state properties, the information and correlations of the 
initial conditions are assumed to be lost. This is archived 
by sending the initial time to — > — oo and dropping all 
correlation functions which involve the initial state. It 
is again assumed, that the system reaches a steady state 
which is translational invariant in time. Therefore, the 
single-particle Green's function does only depend on the 
difference between the two formerly independent times 
of particle creation and annihilation. The Laplace trans- 
form of the time difference then leads to the formulation 
in frequency space for all Green's functions of the steady 
state. 

In the nonequilibrium steady-state formulation two in- 
dependent components of the contour ordered Green's 
function survive which are chosen to be the retarded and 
lesser Green's functions, G r (uj) and G < (uj) respectively. 
The advanced and greater functions are related via 

G a {u) = G r (w) t (23) 
G>(u) = G<(uj) + G r (u;)-G a (u;) . (24) 

The two relevant Green's functions can be expressed 

apmi 
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g<h = |G;H| 2 [2z/ eff H + £<M 

/eff( w ) = /iM r LM + f R (uj)T R {u)) 



(25) 

(26) 
(27) 



where, again, A(w) = A.r(uj) + A^(oj) are the hybridiza- 
tion functions of the leads, F a (u>) their imaginary parts 
(see Eq. ^) and / a (w) = l/{exp /3(w — fj, a ) + 1} are the 
Fermi functions of the corresponding leads. The retarded 
and lesser self-energies, S r and S < respectively, include 
all correlation effects induced by the Coulomb interaction 
U. accounts for the frequency independent Hartree 
energy shift. 



1. Nonequilibrium self-energy 

Diagrammatic expansions in the Coulomb in- 
teraction of the self-energy^ have been investi- 
gated for syste ms in equilibriu nP^ZH as we \\ as j n 
nonequilibrium.EHHIllH13 ti5 | 79 | su | xhe self _ energies can 

be evaluated either non sclf-consistently, where bare 
propagators are used as inner lines, or in terms of 
skeleton diagrams, where fully-dressed propagators are 
taken into account. 

In this study we focus on three different approxima- 
tions for the self-energy: (A) The bare expansion up to 



second order in U, where Hartree-Fock (HF) propagators 
are used as internal lines. The latter are just the non- 
interacting propagators, but with a shifted level position 
E' d = E d + U/2. This approximation is labeled 2 nd U and 
its diagrammatic representation is schematically shown 
in Fig. [lja) and (b). (B) The self-consistent evalua- 
tion of the second-order skeleton diagram of Fig. [l|a) 
and (b). This approximation is called second Born ap- 
proximation (2BA) but in contrast to the usual 2BA, no 
exchange contribution exists for the single impurity An- 
derson model ([2]) with only one spinful orbital. (C) In 
the GW approximation^ 1 ^ (GWA) the bare Coulomb 
interaction U is screened by an infinite scries of particle- 
hole excitations, which can be summed as indicated in 
Fig. [ljc). No contributions with odd orders in the in- 
teraction occur in this series due to the definition of 
the matrix elements of the Coulomb interaction in our 
model ([2]), where we set matrix elements between elec- 
trons with the same spin explicitly to zeroPSl 

The 2BA and the GWA are both evaluated self- 
consistently, and thus the self-energies can be derived 
from a Luttinger-Ward functional!^ Therefore, both 
constitute conserving approximations in the sense of 
Kadanoff and BaymP2 It can be shown that elementary 
sum rules such as charge and current conservation are 
obeyedP3 In contrast, the non self-consistent 2 nd U ap- 
proximation is not conserving which can lead to the vio- 
lation of current conservation, as it will be demonstrated 
later. 

The Hartree shift is produced by the average occupa- 
tion of the quantum dot 



= U{n- a ) 



duj 



2-ni 



(28) 
(29) 



E» = i 



and analytic expressions for the self-energies read 

dx 
2, ( 

dx 

dx 
— ( 
2tt 



+i 



: G<(x)W^u,-x) 
' : G r a {x)W>(cj-x) 
£<M = * J ^G<(x) W<{uj - x) 
where the effective interactions are given by 



(30) 
(31) 



W r (, A - u2p S(^ 

W<(u) =W;(w)P<(w)W»(w) 

and the particle-hole bubbles are 
dx 



(2BA) (32) 

(GWA) (33) 

(34) 
(35) 



2tt 
dx 
2^ 



G r Ax) G<(x 
G<(x) G«{x 



(36) 
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FIG. 1: (a) Schematic diagrammatic representation of the Kadanoff-Baym-Keldysh self-energy. The first term represents the 
frequency independent Hartree shift while the second contribution represents the interaction part. The double-dashed line 
is the effective interaction W while the single-dashed line represents the bare interaction ~ U . The second-order diagram is 
shown in (b). For the non self-consistent second-order U (2 nd U) approximation the internal solid lines with arrows are taken 
as the Hartree-Fock (HF) propagators, while for the second Born approximation (2BA) the internal lines denote fully-dressed 
propagators. In the GW approximation (GWA) (c) the interaction is renormalized by an infinite series of particle-hole pairs 
which can be summed as indicated in the last line. The internal lines again denote fully-dressed propagators. 



, I ~G«{x)G<(x-u) 
, I ~G<{x)G r a {x-u) 
, I ^G<(x)G>(x-u J ) 



(37) 
(38) 



In the above expressions the advanced and greater 
Green's functions can be determined via Eq. (23) and 
p4| and a 



-a denotes the opposite spin of a. 



Equations (25 1- (38) form a closed set, which is solved 
self-consistently for the 2BA and GWA. For the 2 nd U 



approximation all particle-hole propagators (36)-(38) are 



evaluated only once with bare Green's functions 



1 



■iS - E d 



AM 



? <M = 2%»| 2 /effM 



(39) 



(40) 



and Eq. (|32j) is used as the effective interaction. The 
Hartree shift is included in order to determine the desired 
filling. The effective Fermi function f e ff(u>) was defined 
in Eq. ([27]). 

The GWA'SUS2l nas been successfully applied to over- 
come some shortcomings of local-density calculations 
and estimate the screening of the Coulomb interac- 
tions in solid state physics. Recently, it has been em- 



ployed to calcul ate quantum transport through nanoscale 
devices!^ * 20 * 80 ! In the context of the single impurity An- 
derson model it was shown to accurately describe the 
equilibrium properties in the weakly-interacting regime 
and in asymmetric situatio ns w ith a nearly empty or 
nearly full impurity orbital EMS l n the strongly inter- 
acting Kondo regime, i.e. T — U < E d < — V, the GWA 
produces a narrow peak in the spectral function at the 
Fermi level, which could be interpreted as remnants of 
the expected many-body resonance.^! However, the line 
shape of this low-energy resonance as well as the high- 
energy Hubbard peaks at u> sa Ed and u> ss Eg + U a re not 
correctly reproduced by this approximation.^^ Addi- 
tionally, for very large interactions strength U/T > 8, all 
three perturbative approaches favor an unphysical mag- 
netic ground state, which is actually forbidden by the 
Mermin- Wagner theorem.^ In the nonequilibrium situa- 
tion, the proximity to bifurcation points of these sets of 
equations leads to unphysical hysteretic response.^ 



D. Current as function of the bias voltage 



III. RESULTS 



The current flowing from lead a onto the impurity re- 
gion can be expressed as^ 



la = j-S J duT a (u>) \2iG<(u,V) (41) 
+ f a (w)4^p r a {u,V) 



where p r a {u, V) = -3m[G r (u;, V)]/ir is the frequency and 
voltage dependent spectral function of the retarded im- 
purity Green's function. Since the steady-state current 
onto the interacting region from the left must be equal to 
the current leaving to the right lead, i.e. Il = —Ir = I, 
we can symmetrize the left and the right currents with a 
linear combination 66 ! a nd write it as 



I = 



t 2 r Il - rll 



L + R 



(42) 



In the wide band limit, T a {oS) —> T a — r Q (0) and r\ = 
r Q /r holds such that the term proportional to G<(ui) 
drops out of ( 42 ) and we obtain 



I = 



Go 



J2jd" [/lH - /«M] nTpKcj, V) . (43) 



where we have defined Go 



G 



4£lX_r 

r 2 



(44) 



Go reaches the universal conductance quantum e 2 /h for 
a symmetric point-contact junction, — Tr, and is 
strongly suppressed in the tunneling regime T a -C r_ Q . 

For the voltage drop across the two contacts of the 
impurity to the leads we employ a serial resistor model 
where the chemical potentials in the leads are given by 
PL = ~r R V and p R = r\V . 

At zero temperature, the zero bias conductance G = 
edI/dV\v=o = Go7iT P r (0) ^ s proportional to the 
spectral function at the Fermi level. In the zero temper- 
ature Fermi liquid and for a symmetric junction p^(0) = 
l/(7iT). The conductance is given by its universal value 
G = 2Go which shows in the slope at zero bias of the IV 
characteristics, i.e. ie/Go = 2V. 

We also define a leakage current 



AI = 



Il + Ir 
2e 
~h 



4e 

'17 



(45) 



iG<(ui) 



which must vanish due to current conservation, Il = 
—Ir, in a physical junction. Therefore, deviations from 
AI = measures shortcomings of an approximation. 



In this section we compare and discuss the results 
obtained from the different diagrammatic Keldysh ap- 
proaches with the SNRG. For simplicity, we used sym- 
metric structureless leads characterized by a constant 
density of states with a half-bandwidth D = 20T, i.e. 
T a (oj) = T a 0(D - M). The total T = T L + is used 
as the energy scale: All energies, voltages and temper- 
atures are measured in units of T = 1 throughout the 
paper. 

For the SNRG a rather large A — 4 was chosen, and 
N s = 2200 states were retained in each NRG-iteration 
step, z-averaging 6 ^! with either N z = 2 or N z = 4 differ- 
ent z-values was performed, and the broadening param- 
eter for the spectral functiorP^ was chosen b — l.3/N z . 
For large U some unphysical wiggles may emerge in the 
spectral function, as it is explained below. In principle, 
these wiggles can be minimized by choosing a smaller A, 
incorporating more states or performing the z-averaging 
with a larger number of z-values. 

We did not include an external magnetic field, and 
no magnetic solutions are encountered for the parame- 
ter values used in this paper. Therefore, we will drop 
the spin-index from now. The two spin components of 
the spectral functions and self-energies are identical, e.g. 
p r a = p r s = p r and Y7 a = Y7 S = S r respectively. 

Before we apply finite bias voltages, we will compare 
the different equilibrium low-energy scales obtained with 
the diagrammatic approaches to NRG results. While the 
diagrammatic approach becomes exact only in the weak- 
coupling limit U — > 0, the NRG produces the correct 
scales for all interaction strengths. We will identify the 
validity range of the diagrammatic expansion. In that 
regime the diagrammatic approach produces correct re- 
sults even in nonequilibrium, and we will therefore use it 
to benchmark the SNRG for finite voltages. 



A. Equilibrium low-temperature scales 

The single impurity Anderson model in eq uilibr ium for 
T -> always forms a local Fermi liquidPHSEH xhe 
spectral function for a symmetric junction approaches 
the zero temperature limiting value p(u) = 0, T = ) = 
1 / (7rF) in accordance with the Friedel sum rule! 90 * 91 ] The 
Fermi-liquid formation is associated with a characteris- 
tic low-energy scale, which is identified with the Kondo 
temperature Tk at large Coulomb repulsions and near 
half-filling. 

The SNRG coincides with the usual NR d 26 ' 31 ' 88 ' in 
equilibrium, which accurately describes the crossover 
from high to low temperatures and provides the cor- 
rect low-energy scale Tk depending exponentially on f/P^ 
The 2 nd U approximation, however, predicts a low-energy 
scale which is perturbative in U and too largeP^l The 
GWA does produce a narrow many-body resonance in 
the spectral function at the Fermi level. Extracting a 
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low-energy scale from the full width at half maximum 
(FWHM) for an asymmetric junction (Ed ^ —U/2), as 
shown in Fig. 5 of Ref. [551 suggests an exponential varia- 
tion with the ionic level position Ed- However, the expo- 
nent has the wron g prefactor as compared to the exact 
analytic f orm J£2123] 

In order to extract the low-energy scale from our 
model calculations we employ two different methods: We 
calculate the temperature dependent zero bias conduc- 
tanc e G = dl/dV \v=o an d fit it to a phenomenological 
f orm pm Since q j s directly determined by the spectral 
function, it is sensitive to the amount of spectral weight 
in the temperature window —T < uj < T. The scale 
recharge ex t rac t ec j m ^his wa y constitutes the energy scale 
relevant for the zero-bias charge transport in the system. 
This procedure yields the same result as the aforemen- 
tioned extraction from the FWHM of the resonance at 
the Fermi level. 

The second way utilizes the screening of the effective 
local magnetic moment, Me// = Tx(T) = TdM/dH\u—o, 
where x is the magnetic susceptibility, M the magnetiza- 
tion and H an external magnetic field. We calculate M 
for a finite but small external magnetic field SH = 10~ 9 r, 
and extract the susceptibility via the difference quotient. 
In the Fermi-liquid regime the effective magnetic mo- 
ment follows an universal curve as function of temper- 
ature from which the low -energ y scale is determined by 
defining T KX (T K ) « 0.07PMT The resulting T™° s sets 
the scale for magnetic excitations in the system and is 
directly linked to the Kondo- screening of the local mag- 
netic moment. 

For large values of the Coulomb interaction, the scales 
T^ ag and x^ large should coincide (apart from a constant 
of order one) and vary as exp(— nU/8T) for a symmet- 
ric quantum dot. For very small values of the Coulomb 
repulsion U -C T both should approach V. The charge 
scale T'^ lar9e is expected to be roughly constant and on 
the order of T^ har9e ~ T for U/wT < I since for such 
small interactions charge fluctuations to and from the 
leads dominate the physics, and the spectral function 
stays very close to its HF form. On the other hand, the 

m agnet ic scale is known to decrease exponentially for all 
IT mm 

Figure [2] shows the two scales extracted from NRG 
and GWA calculations for a symmetric junction in 
equilibrium. The NRG results show the expected U- 
dependencies: The charge scale T^ arge is on the order of 
r for small U < ST and decreases exponentially for large 
U > 4T. The magnetic scale J 1 ™" 9 decreases exponen- 
tially for all U as it is evident from the comparison with 
a fit function aexp(— ttU/8) also included in the plot. 
Furthermore, there exists only one universality scale for 
large U which manifests itself by T^ ag cx T^ arge (not 
shown) . 

On the other hand, the scales obtained from the GWA 
agree with the NRG only for small U. The charge scale 
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FIG. 2: The equilibrium low-energy scales as functions of U 
extracted from the various approximations as described in 
the text. A fit to the magnetic scale of the NRG to show the 
exponential decay cx exp(— ttU/8) is also included in the plot. 




FIG. 3: Comparison of the NRG and the Keldysh GWA 
equilibrium (V — 0) zero-temperature spectral functions for 
U/T = 4 and Ed = —U/2 and a quantum-point contact 
r_L = r_R = 1/2. The inset shows the temperature evolution 
of the spectral function right at the Fermi level, p r (ui = 0,T). 
The NRG-parameters are A = 2, N s = 1500, N z = 4, 
b = 0.325 and 50 NRG iterations were performed. 



Significant deviations are observed for larger U, where 



the GWA-T^ harge decreases faster than the NRG. For U 
significantly larger than the ones shown in the plot, no 
scales could be extracted due to the artificial symmetry 
breaking already reported in the literature! 80 * ' 

We added a second GWA charge scale j , ^ c ar9e to the 
graph which is obtained from the width of the low-energy 
feature at 75% of p r (0) (and not at the FWHM as for 
rpchargey rpj^ correS p 0nc Ji n giy extracted scale should co- 
incide with Tt ar9e , apart from a prefactor. But it is 
found that both scales follow the same trend only for 
small U and already for U > 3 a much stronger decrease 
than the expected exp(— wU/8T) is observed in r f'^ large . 
Therfore, the extraction of the charge scale within the 



T. 



charge 



K 



perfectly agrees with the NRG curve for U < 4. 
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GWA at intermediate U is somewhat ambiguous. A 
comparison of the zero-temperature equilibrium spectral 
function of the GWA and NRG for U — YY is depicted in 
Fig. [3j The low-energy feature of the GWA spectral func- 
tion is too narrow and exhibits a rather spiky line-shape 
which suggest at too low charge scale. This is supported 
by the evolution of p r (0) as a function of temperature 
which is shown in the inset of Fig. [3j The logarithmic in- 
crease of p r {0) which occurs at temperatures on the order 
of the relevant charge scale also reveals that the charge 
scale is predicted as too low in GWA compared to the 
NRG. However, a considerable broadening occurs away 
from the Fermi level which leads to the same FWHM 
for the GWA as in the NRG and consequently the larger 
recharge emer g e g m thermodynamic quantities like G(T). 

The magnetic scale T^ ag extracted from the GWA ex- 
hibits some peculiar {/-dependence. For small U < Y 
the scale agrees with the NRG. However, it develops a 
minimum at U ~ 2.5r and then increases again for in- 
creasing Ul This clearly indicates a failure of the GWA to 
describe magnetic properties for intermediate and large 
interactions. Since the GWA effective moments Me// 
show universality as functions of the dimensionless tem- 
perature t = T/T K a9 for low temperatures (not shown), 
the increase in T^ acj implies a too strong screening of 
magnetic moments. The effective Coulomb interaction is 
over-screened by W (see Eq. (33)). The electrons remain 
itinerant even at rather large U, and the GWA fails to 
capture the atomic limit. Therefore, the magnetic screen- 
ing scale T^ ag remains large in the GWA and actually 
increases with U. 

The scales extracted from the other two diagrammatic 
approximations all coincide with the GWA for small 
U < r. For larger U the 2 nd U approximation produces 
the same difference as the GWA between the charge and 
the magnetic scale, whereas within the 2BA both scales 
decrease with increasing U , but in a polynomial rather 
than exponential fashion. 

We have established that the diagrammatic approach 
produces reliable results for interactions up to the order 
of the hybridization strength U < T which we will use 
in the following section to benchmark the SNRG in that 
regime. 



B. Weak correlation regime: U < Y 

We study the nonequilibrium properties of a symmet- 
ric and an asymmetric junction in the weakly corre- 
lated regime U/Y < 1. We use a very low temperature 
T = 0.006r, which is sufficiently small compared to all 
other scales in the problem so it can be considered as 
T = with impunity. 

For such small interactions the diagrammatic ap- 
proaches and the SNRG yield identical results for all 
voltages. Figure Qa) shows the nonequilibrium spec- 
tral function of a symmetric junction in the quantum- 
point contact regime, i.e. U = —2Ed = Y = 1 and 
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FIG. 4: (a) Spectral function of the retarded Green's function 
and negative imaginary part of the retarded self-energy (inset) 
of a symmetric junction with U = —2Ed = 1 and Yl = Yr = 
0.5 for different voltages. Results for the spectral functions 
are shown for HF, GWA and SNRG, while the self-energy is 
shown for the GWA only, (b) Current as a function of voltage 
for a symmetric (U = —2Ed = 1 and Yl = Yr — 0.5) and an 
asymmetric junction ( U = 1, Ed = —0.25 and Yl = 4r_R = 
0.8). / is normalized to Go/e = h/e (symmetric) and Go/e = 
0.64/i/e (asymmetric) and measured in units of Y = 1. A 
small temperature of T = 0.006 was used for all calculations. 
Parameters for the SNRG calculations are A = 4, N s = 2200, 
N z = 2, b = 0.65 and 12 NRG-iterations were performed. 



Y L = Y R = r/2 = 0.5. For small voltages V/Y < 0.5 the 
spectra are even indistinguishable from the Hartree-Fock 
(HF) result. Only at larger V small deviations around 
the Fermi level as can be observed. 

The imaginary part of the retarded self-energy 
— Qm[T, r (uj)] for that junction obtained with the GWA is 
shown in the inset of Fig. [l|a) for various voltages. The 
overall scale of — 9m[E r (a;)J is much smaller than Y and, 
therefore, the total self-energy E*°* = S r + A is domi- 
nated by the charge-fluctuation scale Y. But the general 
influence of a finite bias voltage can already be observed 
here: The quasiparticle scattering amplitude is increased 
by the interplay between the voltage-induced fluctua- 
tions and the interaction. The characteristic Fermi-liquid 
quadratic minimum in -3m[S r (w)] at the Fermi level 
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is destroyed with increasing voltage and the local Fermi 
liquid prevailing in equilibrium (V = 0) is suppressed 
at large enough bias. The evolution of the minimum in 
-3m[S r (w)] with voltage bears some resemblance with 
a temperature evolution. The quasi-particle coherence is 
destroyed by a finite voltage in a similar fashion as with 
increasing temperature. 

The resulting IV characteristics of a symmetric and an 
asymmetric junction are shown in Fig. [4jb) for U = T. 
The current is normalized to Go /e and measured in units 
of T. The rescaled current always saturates at 27rr for 
large voltages independent of U (not shown), as required 
by Eq. (43). The initial slope at zero voltage of the IV 



0.3 



curve remains unaltered for all values of U , in accordance 
with the Fermi liquid nature of the model at small bias 
and zero temperature. 

While a symmetric junction with symmetric coupling 
to the leads always results in symmetric spectral func- 
tions, p r (u>, V) — p r (uj, — V), and antisymmetric IV char- 
acteristics, I(—V) — —I{V) (see Eq. (43)), an asym- 
metric junction in combination with asymmetric cou- 
pling yields a non-antisymmetric IV characteristics with 
I(-V) ^ -I(V). This is clearly visible in Fig. gb). 
The bias window ranges from p^ = — r f{V to pu = r\V 
and is not symmetric around the Fermi level. In combi- 
nation with the shift of spectral weight to higher ener- 
gies in p r {ijj) — the center of the spectral function is at 
2Ed + U > — this leads to a smaller contribution to 
the current for negative voltages. 

In contrast to the spectral functions, the IV character- 
istics of the SNRG and GWA agree perfectly with the HF 
results for all voltages. The current is rather insensitive 
to the detailed distribution of spectral weight and mea- 
sures only the total amount in the bias window [/xj,, /xr]. 

The SNRG produces the correct results for small val- 
ues of the interaction U < T, and has thus no principal 
limitations. Therefore, the expectation that it is reliable 
at arbitrary interaction strengths as well is warranted. 



C. Intermediate correlation regime: F < U < 10T 

As demonstrated in the previous section the interac- 
tion plays a minor role for small U/T. On the other 
hand, with an odd number of electrons on the quantum 
dot and at very large U/T and r - U < E d < — r, 
the system develops a Kondo effect as T — > (see for 
example Ref. d and S3)) . The SNRG was showrP to 
correctly describe the strongly correlated Kondo regime 
out of equilibrium. The enhancement of the conductance 
in the Coulomb-blockade region was reproduced for small 
bias, and the destruction of the many-body resonance at 
the Fermi level with increasing voltage has been studied. 

In this section we will focus on the intermediate in- 
teraction regime, where correlations become increasingly 
important. Since the diagrammatic Keldysh approxima- 
tions already show deficiencies in equilibrium — see, for 
example, in Sec. Ill A or Ref. |HH — discrepancies will 
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FIG. 5: Spectral functions of the SNRG and the Keldysh 
approaches for a quantum-point contact with Fl = Tr = 0.5 
and U = -2E d = 4 at T = 0.1. The bias voltages are V = 0.5 
(a), V = 2 (b) and V = 7 (c). SNRG-parameters are A = 4, 
N s = 2200, N z = 4, b = 0.325 and 8 NRG-iterations were 
performed. 



extend to finite voltages. 
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1. Nonequilibrium spectral functions of a symmetric 
junction 

The nonequilibrium spectral function for various volt- 
ages and intermediate interaction U = 4 is shown in 
Fig. |5](a)-(c) for a symmetric junction (Ed = —2) with 
symmetric coupling to the leads (r^ = = T/2 = 0.5) 
at T — 0.1. Now all diagrammatic approximations yield 
different results. 

At low voltages, the SNRG and 2 nd U approximation 
reproduce the slight humps at energies u « ±J7 /2 which 
are the first indicators of upper and lower Hubbard satel- 
lites forming at large U. The GWA only produces the 
broad high-energy tails, without the indication of forming 
separate peaks and the 2BA completely fails to produce 
the enlarged spectral weight at high energies. 

As the voltage is raised, the Coulomb interaction 
causes additional dephasing, leading to increasingly 
broadened spectra. However, the 2 nd U approximation 
produces systematically too broad high-energy tails and 
an unphysical plateau around u) = 0, which even devel- 
ops a slight dip as seen in Fig.[5]jc). We attribute this to 
a tendency to overestimate the Coulomb repulsion. This 
might already be guessed from the equilibrium spectral 
functions, where the 2 nd U approximation unexpectedly 
produces the high-energy Hubbard satellites for arbitrary 
large Coulomb repulsion. These are connected to the 
ionic many-body states of the isolated atom which are 
not expected to be described by a second-order pertur- 
bation theory. However, the analytic structure of the re- 



tarded self-energy, Eq. (30), (32), (36) and (39), has two 
direct consequences: (i) For small coupling to the leads or 
large U it favors a S r (w) oc U 2 / (v+i6) behavior. This re- 
sults in a two-peak structure with the peak-positions and 
widths roughly given by ±{7 and T, respectively. (Incor- 
porating a screened and dynamic Coulomb interaction, 
as it is done in the 2BA and GWA, leads to a prefactor 
smaller than U 2 and additional imaginary parts enter in 
the frequency dependence of S r (w). The Hubbard satel- 
lites are then moved to lower energies and broadened.) 
(ii) The Fermi functions entering Eq. (36) through G < 
lead to a narrowing of the integration interval for decreas- 
ing temperature. At zero temperature this always pro- 
duces a vanishing imaginary part of the self-energy at the 
Fermi level, 3mp r (w = 0,T = 0)] = 0,^ given that the 
noninteracting propagators are non-singular at w = 0. 
(This reasoning also holds for the 2BA and GWA). 

The combination of (i) and (ii) gives rise to the two- 
peak structure in the spectral function for large U and 
the emergence of an additional peak at the Fermi level 
for low temperatures which is usually interpreted as the 
Kondo resonance. But in principle there is no justifica- 
tion why the 2 nd U approximation should be reliable for 
large values of U under arbitrary conditions. Already 
for the asymmetric model in equilibrium the phase-space 
argument (ii) does not guarantee the correct description 
of the low-temperature Fermi liquid anymore, and it is 
well-known that the 2 nd U approximation produces un- 



physical resultsP^I Therefore, the large differences to 
all other methods at finite voltages, as it is observed here 
for U/T = 4, is not surprising. 

The SNRG tends to produce additional features in the 
spectral function for large bias voltages at the positions 
of the chemical potentials of the leads, u> ps These 
are the humps visible for larger voltages in the curves 
of Fig. [5j They are artifacts of the NRG discretiza- 
tion and dependent on the broadening procedure of the 
NRG spectral functionsPSl In equilibrium, the NRG only 
resolves spectral information above a cutoff frequency 
> uj c (T) where ui c is on the orde r of the te mperature 
T. The NRG broadening parameter^ 26 ' 68 1 86 1 94 1 are usually 
adjusted such that artefacts are minimized. Additionally 
the spectral function is interpolated between — lu c < lu < 
u) c . This translates itself to the present implementation of 
the SNRG which does not provide spectral information in 
the intervals ^ = 



-LO c 



centered around the 
two chemical potentials. Here, \x + = max{fiL, Pr} and 
fjT = min{/Z£, ^r}. Furthermore, the time-dependent 
NRG introduces additional discretization errors^ which 
increase with increasing value of U . z-averaging over dif- 
ferent discretizations^ improves the spectral functions 
and these artifacts could be removed by adjusting the 
broadening parameter depending on th e voltage. In this 
paper, however, we keep the broadening^SEHHnB] param- 
eter fixed at b = 1.3 /N z independent of the bias and 
performed z-averaging with N z = 2 and 4. 

Let us focus on the different behavior of the spectral 
functions around u = depicted in Fig. [jjja)-(c). The 
height p r (0) is reduced for increasing V, and the spectral 
functions of the SNRG, 2BA and GWA approach each 
other and eventually coincide. For V — 7, the 2BA curve 
is still a little higher compared to SNRG and GWA, but 
at even larger voltages (not shown) it also falls on top 
of the SNRG and GWA. In contrast, the 2 nd U spectral 
function does not approach this large voltage limit, and 
the zero frequency value of the spectral function is consid- 
erably reduced compared to the other approaches. This 
is in accord with an enhanced scattering amplitude at 
the Fermi level, visible in the imaginary part of the self- 
energy depicted in the insets. 

Upon increasing the voltage, the 2 nd U approximation 
does not follow a systematic trend since p r (0) is larger 
than the SNRG at low voltages and smaller at high V. 
The other approximations show systematic deviations as 
the 2BA is always larger than the SNRG, while the GWA 
is always smaller. 

In the present calculations, the temperature T = 0.1 is 
only about on fifth of the equilibrium Kondo temperature 
for these parameter values, i.e. T/Tk ~ 0.2. As discussed 
in Sec. |III A"} the GWA produces a too small charge scale, 
which results in an even higher effective temperature. 
This leads to a reduction of the spectral function at the 
Fermi level in addition to the effect of the small bias. The 
imaginary part of the self-energy is correspondingly too 
large compared to the SNRG, as can be seen in the inset 
of Fig.[5]ja). The 2BA, on the other hand, overestimates 
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FIG. 6: IV characteristics for U = —2Ed = 4 obtained from 
the spectral functions presented in Fig.[5]for a quantum-point 
contact with Tl = Tr = 0.5. The current is normalized to 
Go/e = h/e and measured in units of T = 1. 



the low-energy scale, which explains the trends in p r (0) 
and -9m[E r (0)]. 

Increasing the current through the junction by apply- 
ing a larger bias enhances the charge fluctuations on the 
local orbital. As already mentioned in section |III B| the- 
ses additional fluctuations introduce dephasing^ and de- 
stroy the coherent quasiparticles which constitute the 
low-temperature Fermi liquid. The accompanying de- 
struction of the characteristic quadratic minimum in 
-3m[E r (u)] around lo ~ is observed in the insets. The 
system is driven away from the equilibrium Fermi-liquid 
fixed point, and the spectral functions at the Fermi level 
decreases. At very large voltage \V/Tk\ 3> 1 the coher- 
ent quasiparticles are completely suppressed, as it can 
be seen from the large imaginary part of the self-energy 
around w « 0. The spiky features in the SNRG self- 
energy for V — 7 are due to the aforementioned dis- 
cretization errors and have no physical meaning. 



2. IV characteristics of a symmetric junction 

A Coulomb interaction U/T — 4 leads to a reduction 
of the current compared to its HF value as depicted 
in Fig. [6j This is characteristic for the onset of the 
Coulomb blockade. All approaches predict this reduc- 
tion but slight differences can be noticed. The 2 nd U 
approximation overestimates the Coulomb blockade re- 
sulting in a current which is systematically smaller than 
the SNRG result. Even though the 2 nd U spectral func- 
tion differs strongly from all other approaches for large 
V, this failure to describe the correct single-particle dy- 
namics is concealed in the current as all approaches yield 
identical results. It again shows the insensitivity of the 
current to the detailed distribution of spectral weight in 
p r (uj). The 2BA slightly overestimates the current for 
intermediate voltages, which is again explained by the 
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FIG. 7: Spectral functions for the asymmetric junction, U = 
4, Ed — —1, with asymmetric coupling Tl = 4r.R = 0.8 for 
small (upper) and large (lower) bias, V = ±0.5 and V = ±7, 
respectively. The vertical dashed lines indicate the location 
of the left and right chemical potentials. The inset shows the 
total occupation (n) of the impurity as a function of the bias 
voltage obtained with the diagrammatic approaches. SNRG- 
parameters are the same as for Fig. [5] 

too large low-energy scale T^ lar3e and the accompanying 
underestimation of correlation effects. 

The GWA current merges with the SNRG result for 
V > 2, which — together with the satisfactory spectral 
function for these voltages — indicates a good description 
of the nonequilibrium properties for intermediate to large 
V. 



3. Asymmetric junction 

The spectral function for a quantum dot with a level 
position Ed = — 1, Coulomb interaction U — 4 and asym- 
metric coupling = 4r^ = 0.8 is shown in Fig. [7] The 
asymmetry between positive and negative voltages is di- 
rectly visible in the spectral functions. 

Apart from voltage-induced broadening which was al- 
ready discussed in previous sections, an additional shift 
of spectral weight in p r (uj) is observed with increasing 
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V. While the SNRG moves spectral weight to higher en- 
ergies for positive bias and towards uj = for negative 
V, the diagrammatic approaches produce the opposite 
trends. 

Due to the stronger coupling to the left lead, T L = 
4Tr, the left-moving scattering states dominate the exci- 
tations and the spectral function of the impurity orbital, 
as can be seen from Eq. ^ and ( 10 ). The effective nonin- 
teracting single-particle excitation energy of an a-mover 
is given by Ae a — Ed~p. a , which implies almost symmet- 
ric parameters for the left-movers at the negative voltage 
V = —7 since then Ae^ = —2.4. Therefore, the asym- 
metry of the spectral function is expected to be reduced 
and p r (uj) to be closer to that of a symmetric junction, 
which is indeed observed in the SNRG. 

The diagrammatic approaches underestimate correla- 
tions in the ionic many-body states, and the occupancy 
of the impurity is overestimated for negative voltages. It 
increases almost linearly with negative voltage as can be 
seen from the inset of Fig. [7^b). Therefore, the Hartree 
shift (28) also increases, and spectral weight is moved 



towards higher energies opposite to what would be ex- 
pected from the physical argument presented above. As a 
consequence the spectral functions are strongly attracted 

to fl L . 

The effective single-particle excitation energy of a left- 
moving scattering state for a positive voltage V = 7 is 
greater than zero, Acl = 0.4. This produces an inter- 
mediate valence situation for the left-movers, where cor- 
relations renormalize the effective excitation energies to 
even larger frequencies^ and a shift of spectral weight 
to higher energies results. An additional drag of spec- 
tral weight towards the chemical potential of the weaker 
coupled right lead, (Xr — 0.8V — 5.6, is expected. The 
SNRG produces such a shift as can be seen in Fig.^b). 
The diagrammatic approaches, however, underestimate 
the level-renormalization in the presence of strong va- 
lence fluctuations, a tendency already observable in equi- 
librium (not shown). Additionally, the reduced occu- 
pancy (inset) diminishes the Hartree energy which again 
leads to a shift towards the stronger coupled chemical 
potential /i£ = —1.4. 

Figure [8] shows the IV characteristics of this junction. 
The asymmetry of /(V) is clearly visible when compared 
to the result from the symmetric junction (also included 
in the plot). For V > 0, the rescaled current is very 
close to its values from the symmetric junction and the 
discrepancies between the 2BA, GWA and SNRG follow 
the already discussed characteristics: The 2BA underes- 
timates correlations and yields a slightly too large cur- 
rent. The GWA has the correct distribution of spectral 
weight in the bias window and produces a rather good 
estimate for the current, despite its deficiencies in the 
description of the single-particle spectra. The failure to 
produce the correct shifts of spectral weights in p r (uj) 
causes the current in the diagrammatic approaches to be 
smaller than the SNRG for negative voltages. 

The 2 nd U approximation, however, reveals its non- 
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FIG. 8: (a) IV characteristics of the asymmetric junction, 
Ed = —1, U = 4 and F L = 4F R = 0.8, calculated with 
the spectral functions depicted in Fig. [7] The result for the 
symmetric junction already shown in Fig. [6] is included for 
comparison, (b) The leakage current AI(V) for the Keldysh 
approaches obtained from Eq. ( 45 I . The currents are normal- 



ized to Go/e — 0.64h/e and measured in units of V = 1. 



conserving nature in the violation of current conserva- 
tion for this asymmetric junction. This is illustrated in 
Fig. ^b) displaying the leakage current AI of Eq. (45 1. 
In contrast to the conserving 2BA and GWA methods 
and the SNRG, AI does not vanish for the 2 nd U approx- 



imation! Thus, left and right current of Eq. (42) do not 
have the same magnitude, i.e. Il =/= —Iri and the cur- 
rent calculated from Eq. (43) does not make sense, since 
different linear combinations ciIl — (1 — o)Ir (0 < a < 1) 
yield different results. Therefore, we did not include the 
calculated IV curves in Fig. [8](a). 

Increasing the asymmetry further, i.e. Tl ^> Fr, 
recovers the equilibrium spectral functions of a quan- 
tum dot coupled to a single lead in all approaches (not 
shown). In the SNRG, the backscattering term 0^ ack , 
Eq. |13[ )) is suppressed, and the model approaches an 
equilibrium single-channel problem. In the diagram- 
matic approaches, the noncquilibrium conditions enter 
only through the effective Fermi function f e ^, Eq. (27), 
which approaches its equilibrium value for F^ — > 0. In 
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this regime, all differences in the spectral functions of the 
presented approaches are given by the known discrepan- 
cies already present in equilibrium. 

IV. SUMMARY 

In the recently developed SNRG approach to open 
quantum systems the scattering states of a noninteract- 
ing quantum impurity model are used to construct the 
nonequilibrium Green's functions for the steady state at 
finite bias voltage. We have established the reliability of 
the SNRG by benchmarking it against the diagrammatic 
Kadanoff-Baym-Keldysh approach, which becomes exact 
in the limit U —> 0. It has been shown that the spectra 
and the current-voltage characteristics agree excellently 
for small Coulomb interactions for symmetric and asym- 
metric junctions at arbitrary bias voltage. 

For intermediate values of U we have compared 
the SNRG to three different approximations obtained 
from the Keldysh approach, namely the second-order 
perturbation theory (2 nd U), the fully self-consistent 
second-order (2BA) and the GW approximation (GWA). 
As correlation effects play an increasingly important 
role discrepancies occur between the different methods. 
These were explained by the insufficient treatment of 
the Coulomb interaction within the diagrammatic ap- 
proaches. 

The Fermi liquid at zero bias voltages is character- 
ized by a single low-energy scale which is captured accu- 
rately by the SNRG, but is not properly reproduced by 
the diagrammatic approaches. No single low-energy scale 
can be extracted from the 2 nd U approximation and the 
GWA at intermediate and large U. While the scale asso- 
ciated with charge fluctuations decreases with increasing 
U, the magnetic scale exhibits a qualitatively different 
U -dependency, since it develops a minimum for interme- 
diate U and increases again towards larger Ul The GWA 
shows a tendency to over-screen magnetic moments for 
increasing values of U and fails to reproduce the atomic 
limit. These deficiencies translate themselves to finite 
bias and explain the discrepancies at small to intermedi- 
ate voltages. 

At large bias voltages the self-consistent diagrammatic 
approaches (2BA and GWA) reproduce the SNRG spec- 



tral functions for a symmetric junction, while the second- 
order perturbation theory yields an unphysical plateau 
around the Fermi level. 

All diagrammatic approximations and the SNRG cap- 
ture the onset of the Coulomb blockade in the IV charac- 
teristics of the symmetric junction. The small discrepan- 
cies are explained by the deficiencies in the treatment of 
the interaction. However, the failure of the 2 nd U approx- 
imation to correctly describe the single-particle dynamics 
at large bias is masked in the current, since there only 
the total spectral weight in the bias window enters. 

In contrast to the other methods, the 2 nd U approxima- 
tion reveals its non-conserving nature by producing a fi- 
nite leakage current for an asymmetric junction, which is 
unphysical. This raises the question about the reliability 
of t he results obtai ned within that method or extensions 
f i t | 22 | 23 | 25 | ioo | ioi | even f or a S y mme tric junction. They 

are only well-justified for cases where |£ r (aj)| <C T for all 
frequencies. 

The voltage dependent redistribution of spectral 
weight for an asymmetric junction is not well-reproduced 
by the diagrammatic approaches. This has been at- 
tributed to too large Hartree shifts due to the wrong 
occupation number of the impurity and the inac- 
curate renormalization of the single-particle level in 
intermediate- valence situations. It leads to the under- 
estimation of the current for large negative voltages. 

The SNRG provides access to the description of 
nonequilibrium steady-state properties of nanoscale junc- 
tions for arbitrary Coulomb interaction and voltages. It 
opens promising perspectives for future investigations, 
such as the influence of charge fluctuations when ap- 
proaching the strongly-correlated regime or the effects 
of an applied magnetic field. 
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